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Abstract: We solve second order relativistic hydrodynamics equations for a boost-invariant 1 + 1- 
dimensional expanding fluid with an equation of state taken from lattice calculations of the ther- 
modynamics of strongly coupled quark-gluon plasma. We investigate the dependence of the energy 
density as a function of proper time on the values of the shear viscosity rj, the bulk viscosity (, and 
second order coefficients, confirming that large changes in the values of the latter have negligible 
effects. Varying the shear viscosity between zero and a few times s/4n, with s the entropy density, 
has significant effects, as expected based on other studies. Introducing a nonzero bulk viscosity 
also has significant effects. In fact, if the bulk viscosity peaks near the crossover temperature T c to 
the degree indicated by recent lattice calculations in QCD without quarks, it can make the fluid 
cavitate — falling apart into droplets. It is interesting to see a hydrodynamic calculation predicting 
its own breakdown, via cavitation, at the temperatures where hadronization is thought to occur in 
ultrarelativistic heavy ion collisions. 
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1. Introduction 

In recent years, the comparison between data from experiments at the Relativistic Heavy Ion 
Collider (RHIC) at Brookhaven National Laboratory on the transverse expansion of the matter 
produced in ultrarelativistic nucleus-nucleus collisions with nonzero impact parameters |IJ, and 
calculations done using second-order relativistic viscous hydrodynamics [H, f|, |5|, EL 171 |^, |S], [IDL 



TT], |T_2|, 13, |14| have strengthened the case that the quark-gluon plasma in QCD at temperatures 
above, but not too far above, the crossover from a hadron gas is a strongly coupled liquid. These 
comparisons indicate a shear viscosity to entropy density ratio 77/s that is within a factor of a few 



of l/47r [ITR M, which is the value of this ratio in any gauge theory with a dual gravity description 



in the limit of infinite coupling and infinitely many colors |16], [17|, |18| . The degree of success of 
the hydrodynamic description of these collisions, which turn out to be creating exploding droplets 
of a fluid that is closer to the ideal liquid limit than water is by about two orders of magnitude 
- and water is the liquid that hydrodynamics is named after - - have refocused attention on 
the question of when a hydrodynamic description applies and when it breaks down. The aspect 
of this question that has drawn most attention is "Before what early time in the collision is a 
hydrodynamic description invalid?" This is the question of how, and how quickly, approximate 
local thermal equilibrium is attained. We shall focus, instead, on the complementary question 
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"Assuming that a hydrodynamic description is valid starting at some early time, after what late 
time does this hydrodynamic description break down?" 

Ultimately, heavy ion collisions produce thousands of hadrons flying outwards toward the de- 
tector. The question of how, and when, a hydrodynamic description breaks down and how one then 
ends up with a cloud of hadrons flying apart is the question of how, and when, "freezeout" occurs. 
In calculations, freezeout is typically handled in one of two ways. One option is to choose (by hand, 
or via fitting the output of the calculation to data) a freezeout temperature well below the crossover 
temperature, and at the surface in space-time where the fluid in the hydrodynamic calculation cools 
to this temperature apply the Cooper- Frye prescription [ID] for mapping hydrodynamic volume el- 
ements directly onto a phase space distribution of noninteracting hadrons. The second option is 
to choose a temperature just below the crossover at which to map the hydrodynamic description 
onto a hadronic transport code, and then let this code describe hadron-hadron interactions until 
the hadrons later freeze out [p0| . The second option is an improvement on the first, but it would 
be even better if the hydrodynamic description itself would tell us when it breaks down. This is 
impossible within the framework of ideal (zero viscosity) hydrodynamics. An ideal hydrodynamic 
evolution is fully specified by its initial conditions and by a thermodynamic equation of state p(e) 
giving the pressure in terms of the energy density. Given p(e), an ideal hydrodynamics code will 
blithely let the expanding fluid evolve until e — » 0, without ever giving any hint that in reality it has 
gone far beyond the epoch when hydrodynamics is actually a good approximation. Thus, one must 
either add a freezeout prescription by hand, or choose by hand to switch from a hydrodynamic to a 
transport description of the expanding fluid. Once viscous effects are included in the hydrodynamic 
description, however, it is possible for the hydrodynamic equations to tell us "from within" when 
they break down. 

There are many ways in which a hydrodynamic evolution of an expanding fluid may break 
down, but we shall focus only on one: cavitation. In an ordinary flowing liquid, cavitation is the 
formation of bubbles of vapor in regions where the pressure of the liquid drops below its vapor 



pressure pl| . Cavitation is well studied, both experimentally and theoretically. It can occur on 
the trailing edges of pump and propeller blades; in this context, the challenge is often to design the 
blade so as to avoid cavitation since when the bubbles produced by cavitation later collapse they 
make shock waves that can damage the blades. Cavitation is also used in medicine — one of the 
treatments of kidney stones involves destroying them via cavitation induced by an ultrasonic pulse. 
Returning to our context, we shall ask whether the pressure in the expanding droplet of quark- gluon 
plasma can become negative, which means that it goes below the pressure of the vacuum — our 
analogue of the vapor phase. This is impossible in ideal hydrodynamics: the thermodynamic p(e) 
is positive. But, in an expanding fluid in the presence of shear and bulk viscosity, shear and bulk 
stresses make an additional contribution to the space-space components of the stress energy tensor 
T^ v and these shear and bulk stresses can drive the pressure — which we shall denote P — negative. 
In particular, bulk viscosity is, in essence, a drop in the pressure of an expanding fluid relative to 
the pressure of an equilibrium fluid at the same energy density, and this decrease in the pressure 
will play a key role in our considerations. We shall only analyze boost-invariant 1 + 1-dimensional 
expansion of a 3 + 1-dimensional fluid [22|, meaning that the fluid is only expanding in the z- 
direction. We shall ask whether, and when, the shear and bulk stresses are sufficient to drive the 
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longitudinal pressure (which we shall denote P^) negative. If P^ goes negative in the hydrodynamic 
equations that describe boost-invariant expansion, what will happen when the fluid expands to the 
point where P^ = 0? At this point, instead of continuing to expand, dilute, and cool in a spatially 
uniform boost-invariant fashion as before, the fluid will break apart into fragments, separated from 
each other by vacuum. Vacuum regions are the analogue of the vapor bubbles in conventional 
cavitation, and their pressure is zero. So, when P^ = regions of fluid can stably coexist with 
regions of vacuuum. The fragments of fluid formed via cavitation will then fly apart, separating 
from each other, and will subsequently hadronize. Our boost invariant calculation will only allow 
us to gauge whether and when P^ goes negative: once the fluid cavitates, it is no longer boost 
invariant, and so our calculation will not be able to describe the subsequent dynamics. Describing 
the size distribution of the fragments that results from cavitation would first of all require including 
transverse expansion in the hydrodynamic description, and would second of all require estimating 
the surface tension associated with the interface between fragments of quark-gluon plasma fluid and 
vacuum. If this surface tension is large, large fragments will result; if it is small, smaller ones will 
be favored. 

We find that as long as rj/s is in the vicinity of 1/Att the shear stress alone is not enough to 
trigger cavitation. 1 However, there is now evidence from a variety of directions [^5[ |28|. |29|, [30 



that the bulk viscosity £ is large — (/s ~ 0(1) 3> r]/s — in a narrow range of temperatures around 
T c . We shall quantify how large this peak in (/s must be if it is to trigger cavitation when the 
expanding fluid has cooled to a temperature near T c . We find that as long as this peak is between 
1/4 and 4 times as wide as suggested by current lattice calculations (whose uncertainties we shall 
discuss), cavitation will occur if the peak is higher than a threshold height that lies between 1/2 
and 1/4 that suggested by the lattice calculations. 

Our paper is organized as follows. In Section 2 we set up the hydrodynamic equations that 
describe the boost invariant 1 + 1-dimensional expansion of a 3 + 1-dimensional fluid, working to 



second order in derivatives of the velocity field [31, 15, 32, 33, 34, 35, p6l pi 37, 38]. After setting 



the full problem up, in the remainder of Section 2 we set the bulk viscosity to zero. We describe 
how we specify the equation of state (which arises at zeroth order), shear viscosity (first order), and 
the various coefficients that arise at second order, as well as the initial conditions. We then show 
results and explore their sensitivity to the second order coefficients and to the shear viscosity. We 
find much greater sensitivity to the shear viscosity, indicating that, as other authors have found 
previously, we are using the second order equations in a regime in which the second order effects 
are much smaller than the first order effects. We turn bulk viscosity on in Section 3. We first 
describe how we parametrize (/s and the one new second order coefficient that at a minimum 
must be introduced and in so doing mention some of the uncertainties in our current knowledge 
of both. We then present our results. In Section 4 we speculate about their implications. The 
possibility that bulk viscosity could cause the expanding fluid to break apart into fragments has 
been discussed previously by Torrieri, Tomasik and Mishustin [RU| using the formalism of Ref. ED 



We close by sketching several facets of the observed phenomenology of heavy ion collisions that 
could indicate that freezeout is preceded by cavitation, some previously highlighted in Refs. 



1 Well below the crossover temperature T c , in the hadron gas, rj/s rises significantly |2^, g4|. If the bulk stress 
does not do so earlier, the shear stress could trigger cavitation at some temperature well below T c . 
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and some not, some coming from recent analyses of data and some of long standing. 2 Because we 
are neglecting transverse expansion throughout we will not be able to make quantitative contact 
with data. But, our results motivate the importance of including the peak in the bulk viscosity 
near T c in hydrodynamic calculations that do include transverse expansion, and the importance of 
looking for cavitation in these more realistic calculations. 

2. Second order relativistic hydrodynamics for a boost-invariant 1 + 1- 
dimensional expansion 

2.1 Setup 

The energy momentum tensor for relativistic hydrodynamics can be written as 

= eu"u u -pA» v + W V (2.1) 

where e and p are the fluid energy density and pressure, is the fluid four-velocity, normalized 
such that = 1, IP" is the viscous tensor, satisfying u^Jl^ u = 0, and where the projector 

A"" = - v?u v (2.2) 

is also orthogonal to v?. Hydrodynamics is the effective theory describing the long- wavelength dy- 
namics of the energy density and the fluid velocity. Its evolution equations describe the conservation 
of energy and momentum, and are given by 

where is the geometric covariant derivative. We shall only be interested in hydrodynamics in 
flat spacetime, but it will be convenient to use curvilinear coordinates to describe boost invariant 
expansion and we therefore keep the geometric notation. With T^ v as in ( |2.1| ), the evolution 
equations take the form 

(e + p) Dv? = V^p - A^DpU ^ , 

Ds = - (e + p) VX + -n^V^) , (2.4) 

where D = u^D^ is the comoving time derivative in the fluid rest frame, V M = A^ U D V is the spatial 
derivative in the fluid rest frame, and the brackets (. . .) denote the combination that is symmetric, 
traceless, and orthogonal to the fluid velocity, namely 

A { ,B V) = (a«A? + A«AJ - l** A„) A a Bp . (2.5) 



2 Explorations of possible consequences of bulk viscosity in the phenomenology of heavy ion collisions that are not 
related to cavitation can be found in Refs. |}7], 
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In general, IP" includes the physics of shear viscosity, bulk viscosity and thermal conductivity. 
Thermal conductivity is only relevant if there is a nonzero density of some species with a conserved 
particle number, and we shall work at zero baryon chemical potential throughout. In a confor- 
mal fluid, the bulk viscosity vanishes and, including terms up to second order in gradients, IP 1 ' 
satisfies |36|, E3 



IP" = rjV {tl u v) - t£ 

_ Ai 
2r) 



A^DIl^ + -IP"V a u c 
3 



(2.6) 



where ui^ u = — \ (V^u u — V u u^) is the fluid vorticity, where the shear viscosity r\ is the only property 
of the fluid that enters at first order in gradients, and where Tq, Ai, A2 and A3 are the four properties 
of the fluid that arise at second order. Obtaining a closed set of evolution equations requires 
specifying the equation of state pie) and specifying rj, rn and Ai^ in terms of e. Equivalently, p, e, 
77, tji and Ai 2,3 can all be specified in terms of the temperature T. It is sometimes also convenient 
to introduce the entropy density 

p + e 



T 



Conformality determines the equation of state p 
-V ^ t>-i 



(2.7) 

e/3 and implies that p = | oc T 4 , 77 oc T 3 , 
rjj oc T _1 , and Ai^.3 oc T 2 , but conformality alone does not determine any of the dimensionless 
proportionality constants other than the one in the equation of state. 

If we relax the assumption of conformality (while continuing to assume throughout that there 
is no net baryon density) the only new term that arises on the right-hand side of ( p.6[) that is 
first-order in derivatives is — C(V a M a )A At!/ , where ( is the bulk viscosity. At second order, many new 
terms arise H2J . As we shall discuss below, it is a standard simplifying assumption to write only the 
term +(T^ l DiV a u a )A tJiU , where r£ is a new second order coefficient whose role we discuss below. 

We shall only consider solutions to the 3 + 1-dimensional hydrodynamic equations in which 
no quantity depends on the transverse spatial directions x and y (which in particular means no 
vorticity) and in which the expansion in the z-direction is boost invariant fl2~2"| . This makes it 



convenient to change variables from (t,z) to (r, £) where r = \Jt 2 — z 2 is the proper time and 
£ = ArcTanh (z/t) is the spacetime rapidity. These curvilinear coordinates are comoving with the 
fluid, meaning that u T = 1 and the spatial components of u all vanish. And, in these coordinates 
boost invariance implies significant further simplifications: IP" is diagonal, and therefore so is T M ", 
and the diagonal components of T^ v depend only on r, not on £. Upon making these simplifications, 
the stress energy tensor in (r, x, y, £) coordinates takes the form 0, [15], | 



33, 



r £P iV 



( e 0\ 
p 
p 

\0 p) 



+ 



/o 








n + |$ 








n + |$ 



\ 





n- $/ 



(2- 
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where the trace of U a/3 — namely II — and the traceless part of U a/3 — namely $ — denote the 
non-equilibrium contributions to the pressure coming from the bulk and shear stresses, respectively. 3 
For a fluid at rest, the pressure is isotropic and is given by p, which is related to the energy density 
by the thermodynamic equation of state p(s). As the fluid is expanding, unless it is an ideal fluid 
with II a/3 = its pressure is no longer isotropic — the transverse and longitudinal pressure are 
given by 



P± = p + n + -$ 

p 5 = p + n - $ . 

Furthermore, upon making these simplifications the second order evolution equations are 

e + p + n - $ 



T, 



4?7 
37 ' 

c 



n 



3r 



2r] 2 



' n dr ~ \ 

At first order, II and $ are given by their Navier-Stokes values 
and 

n 



C 



(2.9) 
(2.10) 

(2.11) 
(2.12) 

(2.13) 

(2.14) 
(2.15) 



We see that if we ignore the terms in square brackets in Q2.12|) , then the second order equations ( |2.12| ) 
and Q2.13Q describe $ and II relaxing towards their Navier-Stokes behavior (p. 14 ) and fl2.15|) with 



time constants and t-q, along the lines of the Israel- Stewart approach to second order dissipative 



relativistic hydrodynamics [fl~5"] . If we ignore bulk viscosity, setting II = 0, then (|2.11|) and (|2.12f) , 
including in particular the terms in the square brackets in (|2.12|), follow from conformality 



However, once we turn on bulk viscosity we are breaking conformality, and there can then be further 



second order terms in both ( ^.12[ ) and ( ^.13| ) [13]. These equations are in this sense provisional, 



but it should be noted that the terms in square brackets in ( |2.12| ) become neglible at large r and 
we expect the same to be the case for the missing nonconformal terms also. 

2.2 Signs of cavitation 

Since $ > and II < at first order, see (|2.14j ) and (|2.15|) , it is reasonable to expect them to have 
these signs in solutions to the second order equations also. We then see from ( |2.9p and ( |2.10| ) that 
if either ( or rj is large enough, the longitudinal pressure can be driven negative, and if ( is large 
enough, the transverse pressure P± can also be driven negative. We shall see in Section 3 that if 
( rises high enough at temperatures in the vicinity of the crossover from quark-gluon plasma to 
hadron gas, the resulting bulk stress II does drive P^ negative, indicating cavitation. 

3 The quantity that we denote as <J> has been called $ in some of the literature and II elsewhere in the literature. 
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2.3 Equation of state, shear viscosity, and Ai 



We see that in order for the evolution equations (|2.11|) , (|2.12|) and fl2.13| ) to be closed we need the 
equation of state p(e) and expressions relating rj, (, r^, and Ai to e. For the remainder of Section 
2, we shall set ( = and therefore II = 0, deferring our analysis of the effects of bulk viscosity to 
Section 3. We then need p, r], and Ai only. 

We need an equation of state p(e) that describes the quark-gluon plasma phase as well as the 
crossover to a hadron gas. Lattice quantum field theory is well-suited to the calculation of thermo- 
dynamic quantities at zero baryon chemical potential, and so there are many lattice calculations of 



p(e) in QCD that we could employ. We shall take p(e) from Ref. |46], both because it is an example 
of the state of the art and because these authors have provided a parametrization of their results 
that is easy to use. They parametrize their results for the trace anomaly using the functional form 



e — 3p 



1 



\ 



1 + exp 



-ci 



C'2 



2 

J / 



d 2 G?4 
7^2 J*! 



(2.16) 



and give values with error bars for the coefficients d,2, d^, c\ and c 2 for calculations done with two 
different lattice actions, with and without combining these calculations with hadron resonance gas 
calculations valid at lower temperatures. We shall use the central values of their results obtained 
from combining lattice calculations done with the p4 action and hadron resonance gas calculations: 
d 2 = 0.24 GeV 2 , d 4 = 0.0054 GeV 4 , c x = 0.2073 GeV, and c 2 = 0.0172 GeV. These authors 
find a crossover between hadron gas and quark-gluon plasma occurring in a temperature regime 
180 MeV< T < 200 MeV. In Section 3 when we need to specify a value of the crossover temperature 
T c we shall use T c = 190 MeV, in order to be consistent with the equation of state that we employ 
throughout. The pressure is related to the trace anomaly by 



P(T) _ p(T 

y4 



1 



dT' 



3p 



To 



JV5 



(2.17) 



and the results of Ref. jKJ are obtained by choosing T = 50 MeV and p(T ) = 0. We shall only 
work at T > 100 MeV, where there is no effect of these choices. Knowing (e — 3p) and p as functions 
of T, we know e as a function of T also, as well as the entropy density (|2.7|) . And, from p(T) and 
e(T) we have the equation of state p(e). We shall use the same equation of state throughout this 
paper, focussing on the effects of varying other quantities. 
Next, we turn to the shear viscosity 77. We shall use 



V 

s 



1 

47T 



(2.18) 



as a baseline value, and we shall explore the effects of varying rj/s. The relationship ( |2.18| ) holds 
for the plasma phase of any gauge theory that has a dual gravity description, in the limit of large 
numbers of colors and infinitely strong coupling [ITBL ITTTL IT3. Even though much larger values of 



rj/s (of order 1 and larger) are expected both in the hadron gas found well below T c [23, 24] and 
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in the weakly coupled quark-gluon plasma found far above T c j^J, the baseline ( [2.18| ) is seen as 
a reasonable starting point for the analysis of quark-gluon plasma in the regime being explored 
by RHIC collisions, say around T ~ 1.5T C . Lattice QCD calculations, to date in a gluon plasma 
without quarks, indicate values of 77/s only a few times (|2.18|) [ 481 |30| : at T = 1.58T C Meyer finds 
(77 + j()/s = 0.20 with a statistical error of ±.03 PP| . (At this temperature, ( is small compared to 



77.) Comparison between second order relativistic viscous hydrodynamic calculations that include 
transverse expansion H [4], ||, ||, |7|, ||, |9|. [10|, [DJ [L3], [L4] and data from RHIC |L], |2[ on the azimuthal 
anisotropy of ultrarelativistic heavy ion collisions that have a significant impact parameter indicate 
that in these collisions approximate local thermal equilibrium is attained rapidly and 77/s is small, 
apparently < 0.2 and conservatively < 0.5 |15|, [14[]. Given that we analyze longitudinal expansion 
only, we can have nothing to say about the extraction of information about 77/s from these data. 
But, we shall confirm that 1 + 1-dimensional boost invariant expansion is modified significantly as 
we vary 77/s between and 2/4ir. 

Less is known about the values of the second-order coefficients and Ai. We shall take as a 
baseline 

v 2-log2 



2ttT 



(2.19) 



and 



Ai 



rj 



2ttT 8n 2 T 



which are their values in the plasma of M = 4 supersymmetric Yang-Mills theory fl36 



(2.20) 
the 

simplest and best studied example of a strongly coupled plasma with a dual gravity description. 
Meyer finds 2ttT = 3.1 ± 0.3 at T = 1.58T C in lattice calculations of QCD without quarks [j30] . 
within a factor of a few of (|2.19|) . 4 We shall find that the effects of varying and Ai by large factors 
are small, indicating that the hydrodynamic calculations are being done in a regime in which second 
order effects are small compared to first order effects. 



2.4 Baseline results 



In Fig. [I] we show an example of a solution to the evolution equations ( |2.11| ) and (|2.12| ) with 
vanishing bulk viscosity. The pressure p and temperature T are related to the energy density e via 



the lattice calculations of QCD thermodynamics from Ref. [M\ that we have described in Section 
2.3. The shear viscosity 77 and the second-order coefficients and Ai have been set to their baseline 
values (|2TT8p , Q2Tg ) and (gj§ . 

We have initialized the evolution in Fig. [TJat a time r = 0.5 fm/c, a reasonable choice given that 
the RHIC data on the anisotropic expansion in heavy ion collisions at nonzero impact parameter 
can only be understood if a hydrodynamic description is already relevant earlier than 1 fm/c after 
the collision p2"| . We have also chosen a value for the initial energy density that is reasonable for 
collisions at the top RHIC energy. At time r = 1 fm/c in the evolution of Fig. [I], the energy density 



4 In kinetic theory, 7y} = (5.0 to 5.9)?7/(Ts) [[30| |3^, pl|| , with the prefactor depending on the value of the coupling 
constant. Kinetic theory is not quantitatively valid if rj/s = 1/47T, but applying it anyway gives TyJ in agreement 



with Meyer's lattice result. In kineti c the ory, Ai = (4.1 to 5.2)ri 2 /(Ts) |5C 
value of Ai within a factor of two of ( 2.20 ). 



which with rj/s = 1/Att would give a 
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320 




x [fm/c] x [fm/c] 



Figure 1: Evolution of the energy density e (plotted as e/3), pressure p, and shear stress 3> (left panel) 
and temperature T (right panel) as functions of proper time r. The equation of state is taken from lattice 
calculations of QCD thermodynamics as described in Section 2.3 and the shear viscosity r\ and and Ai 
have all been set to the baseline values described in Section 2.3. The bulk viscosity is zero. The evolution 
starts at t = 0.5 fm/c with the energy density e being that at T = 305 MeV. The evolution starts with 
$ = 0. 



is e = 7.12 GeV/fm 3 , which is consistent with estimates of the energy density at this time that have 
been made using data on the final state energy and entropy [53|]. (Using the lattice results for 
s(T), this energy density corresponds to a temperature T = 252 MeV.) So, although it is misplaced 
precision to specify initial conditions with e = 15.8 GeV/fm 3 (and T = 305 MeV) at r = 0.5 fm/c, 
we shall make this choice throughout this paper as in Fig. p] since varying these choices within a 
reasonable range would have no qualitative effects. Note also that if we were doing phenomenology 
(which we cannot do given our 1 + 1-dimensional expansion) we would want to adjust the initial 
state energy density as we vary other parameters (which we will do below) in order to maintain the 
same late time energy density. We shall not do this, since our purpose is to explore how solutions 
to the evolution equations depend on parameters, and tweaking the initial conditions as we varied 
the parameters would for this purpose be a complication. 

In Fig. H we have chosen $ = at r = 0.5 fm/c. There is no phenomenological justification for 
this choice. Instead, we find that this choice does not matter. What we observe in the evolution is 
that $ rapidly (over a timescale that is a few tenths of a fm/c in Fig. [I] and that is controlled by 
Tq) increases until it is close to its Navier-Stokes behavior (|2.14f) , and then follows ( p.!4| ) closely 
during the subsequent evolution. If instead of initializing with $ = we choose $ at r = 0.5 fm/c 
to be twice its Navier-Stokes value, we find the same behavior. And, varying the initial value of $ 
over this range makes very little difference — it changes e(r) by less than half as much as we shall 
find when we vary in the next section. 

We have plotted Fig. | up to a proper time of r = 10 fm/c, when T = 171 MeV, well below 
T c ~ 190 MeV. Extending the calculations to later times, we find T = 156 MeV at r = 20 fm/c, 
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4 6 
x [fm/c] 




x [fm/c] 



Figure 2: Effects of varying and Ax- The left panel shows the effects on the energy density e(r) of 
increasing t-q. From top to bottom, the curves show e(r) with increased by factors of 2, 4, 8 and 16 
relative to its baseline fl2,19| ) divided by e(r) with r^J unmodified. The right panel shows the effects on the 
energy density e(r) of increasing Ai relative to its baseline ( [2.20D by factors of (top to bottom) 2, 4, 8 and 
16. 



but this is not relevant because at these low temperatures, the shear viscosity of the hadron gas is 



much greater than the baseline value ( |2.18| ). 

Note that in Fig. [I], the shear stress $ is less than the isotropic pressure p at all times, meaning 
that the longitudinal pressure of (|2.10|) is everywhere positive. We shall see in Section 2.6 that 
this need no longer be so if larger values of i]/s are used. 



2.5 Insensitivity to and Ai 

In Fig. || we see that both the second order coefficients and Ai can be increased by large factors 
without having significant effects on the evolution of e{r). Increasing extends the initial period 
of time in the evolution when the shear stress $ changes from its initial value to approach its 
Navier-Stokes behavior ( [2.14 ). With increased relative to its baseline value ( |2.19| ) by a factor of 
8, these early-time transients last 1 — 2 fm/c. Even in this case, Fig. |2| shows that the modification 
of e(r) is small. Increasing Ai depresses <3> relative to its Navier-Stokes behavior ( |2.14j ) at all times, 
but the effect is small even when Ai is increased relative to its baseline value ( |2.20| ) by a factor of 
16, and the effect on e{r) is even smaller. To see significant effects on $, Ai must be increased by 
factors of ~ 100, and even then the effect on e{r) is only of order 10 percent. We conclude that 
we are using the evolution equations in a regime in which the effects of the second order terms are 
small — smaller, we shall see, than the effects of the first order terms. We shall therefore use the 
baseline values Q2.19|) and (|2.20|) exclusively in results that we shall quote throughout the following. 
But, we have checked that varying these parameters simultaneously with the variations of rj and ( 
that we discuss below does not change any interesting conclusions. 
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Figure 3: Effects of varying the shear viscosity rj. The left panel shows the effects on the energy density 
e(r) of setting rj/s to 10~ 4 (lower curve) or 2/4-7T (upper curve), relative to its baseline value of 1/4tt. (The 
curves show the ratio of e(r) with the modified rj/s to e(r) with rj/s = l/4ir.) The right panel shows e/3, 
p, and the shear stress as functions of r for the case where r]/s = 2/4tt. (This panel should be compared 
to the left panel of Fig. |l[) We see that the longitudinal pressure P^, which in the absence of bulk viscosity 
is given by p — <£, comes close to vanishing at an early time. 



2.6 Sensitivity to shear viscosity 

In Fig. |] we see that changing the coefficient of the first order term in the evolution equations, 
namely the shear viscosity i], has much more significant effects than those we found in Section 2.5. 
Increasing r\ by a factor of two has a ~ 15% effect on the energy density sir). The fact that the 
evolution equations are much more sensitive to variation of 7] than they are to variations of or Ai 
provides qualitative support to the program f|, ||, ||, |7|, [5], |9], [TD|, [TJ], [12], [13], [L4[] of using comparison 
between hydrodynamic calculations that include anisotropic transverse expansion and data from 
RHIC to extract information about the value of rj/s — this extraction should not be complicated 
by our lack of knowledge of the values of r n or A. We will revisit this conclusion in Section 4 after 
considering the effects of bulk viscosity, which we have so far neglected, in Section 3. 

Furthermore, we see from the right panel of Fig. [I] that, as Martinez and Strickland have 
analyzed in detail p8| , increasing r] makes the longitudinal pressure of ( |2.10| ) negative at early 
times. With the initial conditions that we are using, we find that is negative for some window 
of early r's if rj/s > 2.15/47T. Martinez and Strickland have analyzed how this criterion depends 
on the choice of initial conditions. It is easy to see why < must arise at early times for 
sufficiently large rj/s. We have seen that the shear stress rapidly rises to its Navier-Stokes value 
( Enp . Together with rj oc s this means that, after transient behavior at very early times, $ oc s/t. 
The zeroth order (ideal hydrodynamics; no shear viscosity) solution to the evolution equations for 
boost-invariant expansion with an equation of state p oc e has p oc e oc r -4 / 3 and s oc 1/r, meaning 
that if $ oc s/t then $ grows faster for r — > than p does. This means that at some early time, $ 
must exceed p making P^ < 0. Our results and the results of Ref. [38] confirm that the conclusions 
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of this simple argument apply. With the specific initial conditions that we have used — namely 
with $ = initially — negativity of P% can be avoided at any given rj/s by increasing by a large 
enough factor. This stretches out the initial transient, delaying the r at which $ oc s/r is reached 
until late enough that $ never exceeds p. However, this resolution is specific to our (completely 
arbitrary) choice of the initial value of $ — we could have chosen $ = 4^/3r from the beginning 
- and is therefore not actually a resolution. The conclusion we should draw is simply that the 
hydrodynamic description, premised on local thermal equilibrium, must break down at early times 
and the negativity of P^ is one sign that tells us before when we cannot go. (Other evidence for 
the same qualitative and, essentially, quantitative conclusion has been developed in Refs. |35|, [5TJ . 



For rj/s = 1/47T, initializing the hydrodynamic evolution equations at tq = 0.5 fm/c as we are 
doing is appropriate. For rj/s = 2/Air, choosing To = 0.5 fm/c is only appropriate for certain initial 
conditions (including $ = 0), while choosing tq = 1 fm/c is safe for more generic initial values of 
$. For rj/s = 3/47T, we find that in order to avoid P^ < we must choose tq > (3 — 4) fm/c. It 
would be interesting to pursue this line of reasoning further in the higher-dimensional hydrodynamic 
calculations of Refs. [|, [5], f| 0, § || OH 0> HH HH : ^he same hydrodynamic calculations that 
yield information about the allowed values of rj/s via comparison to data from RHIC should at 
the same time constrain the earliest time at which a hydrodynamic description can be valid, via 
checking before what time these hydrodynamic evolutions feature negative pressure in some region 
of space. 

The shear viscosity to entropy ratio rj/s becomes large at late times, in the hadron gas phase ^3 



24[ . In our calculations with rj/s = I /Air and rj/s = 2/4-7T, the shear stress $ is much less than 
p at late times. But, if we consider rj/s increasing at late times to rj/s ~ 1, as appropriate at 
T ~ 100 — 150 MeV according to the calculations of Demir and Bass |f24}l , we find Pg coming close 
to going negative at the very late times corresponding to T ~ 150 MeV. This would be worth 
further investigation, as a possible indicator of when the hydrodynamic description breaks down at 
late times, if not for the fact that freezeout in heavy ion collisions is expected to occur earlier than 
this. And, furthermore, we shall see that including the effects of bulk viscosity can result in the 
breakdown of hydrodynamics also happening earlier, when T ~ T c . To this we now turn. 



3. Effects of bulk viscosity 

We now wish to turn on a nonzero bulk viscosity ( and study its effects on solutions of the evolution 
equations ( |2.11| ), ( |2.12| ) and, now, fl2.13| ). We shall set r/, and Ai to their baseline values ( |2.18| ), 



fl2.19| ) and ( |2.20|) throughout this Section. In the case of and Ai we do so because, as in Section 
2.5, we do not expect that changing their values would have significant consequences. But, we have 
seen that varying rj/s is consequential. We are setting rj/s = 1/Att in order to be conservative, in 



the following sense. We shall focus on the question of whether the longitudinal pressure P^ of ( [2.10 ) 



goes negative. Increasing r] increases which makes a negative contribution to P^. So, if we find 
that the bulk viscosity drives negative with rj/s set to 1/47T, increasing rj/s would only make P^ 
even more negative. 

Both at low temperatures T <^ T c in the hadron gas p4| , |55L |SB|, |57] and at very high temper- 



atures T>T C where the quark-gluon plasma is weakly coupled |58| the bulk viscosity ( is much 
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smaller than the shear viscosity. These calculations indicate that (/s rises as one approaches T c from 
both below and above. And, if the crossover at T c were a second order phase transition, ( would 

The general expectation that (/s may be significant near T c is supported 



29 



peak at T c [^, | 

by the analysis of Refs. |26, ^8f which uses sum rules to relate the bulk viscosity to (derivatives of) 
thermodynamic quantities calculated on the lattice (although this relation is subtle |59], |6(], |61|) and 



by the analyses of Ref. [£7j] i n which ( itself is constrained via lattice calculations, albeit in QCD 
without quarks. Both these approaches find C/ s ~ 1 ^ i]/s in a narrow range of temperatures very 
near T c . We shall ask whether such a peak in (/s can drive negative, triggering cavitation. 



3.1 Choosing (/s and 

Determining (/s via lattice calculations of Euclidean correlation functions is challenging, and the 
results obtained in Ref. |27| should not yet be seen as definitive ]55|, [52], |B3| . To give just one example 
of a difficulty [^|, just as ( peaks at a second order phase transition, so does the relaxation time r£ 



- due to the phenomenon of critical slowing down. The Euclidean lattice calculations are sensitive 
to the ratio C/ r n> making it hard to disentangle one from the other. Note, however, that attributing 
a peak in C/ r n to C is conservative in the sense that by neglecting the rise in one underestimates 
the rise in (. Unfortunately, other more technical challenges can go in the other direction [52, |63 |. 
Progress is nevertheless possible 164], |30||. The current state of affairs is that lattice calculations 



make a robust case for the existence of a peak in (/s at T c , at least in QCD without quarks [[30 



but they do not yet reliably determine the height of the peak. We shall therefore parametrize the 
results of Ref. P7| , and vary the parameters over a considerable range. 

In Ref. |2"7|] , Meyer reports results for (/s at five values of the temperature, T/T c = 1.02, 1.24, 
1.65, 2.22 and 3.22. One way to parametrize his results is to write 



c 



- = aexp 
s 



T c -T 
AT 



for T>T r 



(3.1) 



with a, AT and b the parameters. (T c is not a parameter in that QCD without quarks has a 
first order phase transition at a reliably determined T c . When we employ (|3.1|) together with the 
equation of state for QCD with quarks specified via (|2.16|) , we shall set T c = 190 MeV.) Meyer's 
results at the three higher temperatures, well above T c , are consistent with (/s oc 1/T 2 . This is not 
surprising since the trace anomaly (e — 3p)/(e + p), which like (/s is a dimensionless measure of 
the breaking of conformal invariance, is oc 1/T 2 at high temperature. (For example, see fl2.16p .) If 
we set a = and fit b to Meyer's central values of (/s at the three higher temperatures, we find 



0.061 



(3.2) 



Note that at these higher temperatures, Meyer gives the central values for (/s that we have used 
but his results remain consistent with ( = 0. When we vary b, therefore, we should consider values 
as small as zero. With b chosen as in ( |3.2| ) and with a = 0, the curve ( |3.1|) is far below Meyer's 
results at T = 1.02T C and T = 1.24T C . That is, there is no way to use simply C/s oc 1/T 2 to fit 
Meyer's high temperature results and his results at T = 1.02T C and T = 1.24T C — where (/s has 
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its peak. The parameters a and AT can then be chosen such that ( |3.1| ) passes directly through 
Meyer's central values of (/s at these two temperatures, yielding 

a = 0.901 and AT = -p- . (3.3) 

14.5 

The parameter a controls the height of the peak in (/s and AT controls its width, and we shall 
vary both considerably. 

Little is known about the value of Tq. We shall use 

4 = rfj (3.4) 

with given in ( |2.19| ) as the baseline value for r^. Although there is no strong argument for 
this choice, it holds in one class of strongly coupled nonconformal fluids |6^, [44|. 5 On general 
grounds ( 1.6. clS £L manifestation of critical slowing down) and as in the specific strongly coupled 
nonconformal fluid studied in Ref. [36|, t-q is expected to peak where ( peaks, and we shall therefore 
check the effects of greater than in (^.4[) by as much as a factor of 40. 



3.2 Boost invariant expansion with bulk viscosity 

Once one has picked (/s and r^, the next step is to choose initial conditions. We shall initialize 
at To = 0.5 fm/c and choose e as in Section 2. We shall choose the initial shear stress $ = as in 
Section 2, and we shall also choose the bulk stress II = 0. When we then evolve the equations of 
motion (|2.11|) , fl2.12|) and ( |2.13|) , we find that II quickly evolves to its Navier-Stokes value ( |2.15|) 



during an initial time that is controlled by t£. We shall always stop the evolution at the time when 
T has dropped to T c , since our parametrization of (/s in ( |3.1| ) is only valid for T > T c . We expect 
C/s to drop rapidly below T c , but even less is known about the shape of the (/s peak below T c 
than above it, so we simply stop the evolution when T = T c and ask whether by that time the 
longitudinal pressure P^ has gone negative. 

In describing our results, we begin with a = 0. That is, we begin with no peak in (/s, just with 
(/s oc 1/T 2 . If we choose b as in (|3.2|) , suitable to describe the lattice results at T > 1.5T C |p7 



we 



find that the bulk viscosity has negligible effects. Introducing the bulk viscosity changes the energy 
density by about 6%. And, with b as in ( |3.2|) , we find that we must increase the shear viscosity 
r]/s to 1.8/47T in order to see P^ < at early times — whereas we saw in Section 2.6 that with no 
bulk viscosity this required rj/s = 2.15/47T. If we reduce b relative to (|3~2|), the effects of the bulk 
viscosity become even more negligible. Even if we increase b by a factor of 2 relative to (|3.2j ), we 
still find no qualitative consequences: the energy density changes by about 13% and P^ < at early 
times requires 7]/s > 1.4/47T. Note that with b greater than (|3.2| ) by a factor of two there is already a 
range of temperatures near T c where ( > rj. Increasing b by another factor of two makes ( > rj over 
a wide range of temperatures, which is not supported by the lattice calculations. Henceforth, we fix 
b as in ( |3.2| ), meaning that if a were zero the bulk viscosity would not have interesting consequences. 

We now investigate the consequences of the peak in (/s near T c . Let us begin by choosing a 
and AT as in fl3.3| ), meaning that the peak in (/s has a height and width as indicated by Meyer's 



3 In kinetic theory, = 
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Figure 4: Consequences of a peak in the bulk viscosity near T c . We plot e/3, p, the shear stress <3?, the 
bulk stress IT, and the longitudinal pressure P^ = p + U — $ as functions of proper time r. In all the panels, 
the curves end at the r when T = T c = 190 MeV. But, in all the panels, the peak in the bulk viscosity 
near T c has driven negative at some earlier time. After these times, the curves are not relevant because 
when P^ reaches zero instead of continuing to expand smoothly the fluid would cavitate, falling apart into 
regions of fluid separated by regions of vacuum. In the top-left panel, the bulk viscosity is as in (|3.1|) with 
(3.2) and ( |3.3| ). In the top-right panel, a has been reduced relative to ( |Q| ) by a factor of two — the peak 
in the bulk viscosity near T c is half as high. In the bottom-left panel, a is as in (^^) and AT has been 
reduced relative to (|3.3j) by a factor of four — the peak in the bulk viscosity is four times as narrow. In the 



bottom-right panel, a and AT are as in ( |3.3|) but = 20tjj, whereas in all other panels 7"n = r n- I n an 
the panels, the parameters rj, and Ai are set to their baseline values, as in Fig. la. (In the top-left panel, 
P^ = and cavitation occurs at r = 2.3 fm/c when T = 211 MeV; in the top-right panel, at r 
when T = 195 MeV; in the bottom-left panel, at r 
panel, at r = 3.5 fm/c when T = 197 MeV.) 



4.2 fm/c when T 



3.7 fm/c 
193 MeV; in the bottom-right 



lattice results [p7| . We illustrate the resulting evolution in the top- left panel of Fig. |. We see that 
with a and AT as in ( p.3| ), the rising bulk viscosity drives the longitudinal pressure negative 
when T is still well above T c . As we have described in Section 3.1, although there is good evidence 
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for a peak in C/s near T c , its height and width (which we are parametrizing by a and AT) are not 
well known. So, we have explored for what values of these parameters P^ is driven negative near 
T c . As the top-right panel of Fig. f| shows, upon reducing a by a factor of 0.5 while keeping AT 
fixed we continue to find P^ < 0. In fact, we find that with AT as in (|3.3| ), the largest a at which 
P^ remains positive for all T > T c , which we shall denote a sta bie, is citable = 0.37, which is 41% of 



the value (|3.3| ) indicated by the lattice calculations of Ref. ||27||. How does changing the width AT 
modify this result? If we increase AT by a factor of two relative to ( |3.3| ), a sta bie drops only slightly, 
to 0.33. If we increase it by another factor of two, making the peak in (/s four times wider than in 



Ref. J27|, ^stable drops to 0.22. If we change AT in the other direction, decreasing it by a factor of 



two relative to ( |3.3|) , meaning that we make the peak in £/s twice as narrow as in Ref. |27j, a s tabi. 



rises only very slightly, to 0.39. If we decrease AT by another factor of two, arable rises to 0.45; if 
we decrease it by yet another factor of two — making the peak in (/s eight times narrower than 



in Ref. [£7| ~ Stable increases to 0.54. In the bottom- left panel of Fig. f| we illustrate the case 
with AT reduced relative to ( |3.3| ) by a factor of four, and a unmodified, as in (|3.3j). Comparing 
this panel to the top-left panel, we see that reducing the width of the peak in the bulk viscosity 
by a factor of four delays the time at which P^ goes negative, but does not significantly reduce the 
amount by which it goes negative. This is consistent with the observation that a sta bie is not much 
changed. We conclude that a sta bie is quite insensitive to AT over a wide range of AT, ranging from 
0.22 if AT is four times larger than (|3.3| ) to 0.45 if AT is four times smaller than ( |3.3|) . Over this 



wide range, a sta bi c is well below the value 0.90 indicated by the lattice calculations of Ref. p7 



Note that all the results we have just quoted were obtained with r]/s = 1/4tt. If we increase r], 
it will take a smaller (/s to push the longitudinal pressure negative. For example, with AT as in 



( p.3|) if we increase r]/s to 2/47T this decreases the value of a sta bi e from 0.37 to 0.32. 

It is also worth checking that the fact that P^ is being driven negative really is due to the peak 
in the bulk viscosity, not to the 1/T 2 component of (/s in ( |3.1| ). To this end, we set b = (with 
rj/s = 1/47T and AT as in ( |3.3|) ) and found that eliminating the 1/T 2 component of £/s increased 
Q stabie, but only from 0.37 to 0.41. 

The curves in all the panels in Fig. |] end at the r when T = T c in the calculation, although they 
have ceased to be relevant earlier when P^ reaches zero and cavitation occurs. Note, however, that 
the r at which T = T c is between 4.7 fm/c and 5.2 fm/c in all four panels, while in Fig. [I] T reaches 



T c at 4.3 fm/c. This is a small effect, but it can be made larger. As Kapusta discovered in Ref. |)7 



if (/s diverges at T c , namely if C/s oc 1/|T — T c \ n for some positive power n, then the diverging bulk 
viscosity acts in the hydrodynamic equations so as to prevent T from dropping below T c : as r — > oo 
the equations show T approaching T c from above but never reaching it. The slight delay in the r 
at which T c is reached in our calculations, as in Fig. f|, is the small residue of this effect when the 
peak in (/s is finite. Note that the solution with diverging (/s and T approaching T c from above 



We have used the equation of state (2.16) and (2.17) obtained from lattice QCD calculations throughout. If 



instead we use the conformal equation of state p = e/3 with no phase transition, then a s t a bio = 0.95. It is easy to 
see why a larger peak in the bulk viscosity is then required in order to drive the longitudinal pressure negative: in 



the vicinity of the peak in the bulk viscosity, p = e/3 is significantly larger than p in ( 2.17 ) because ( 2.17 ) describes 
a phase transition; because p is larger, it takes more bulk stress to drive negative. One measure of the robustness 
of our results is that even with the much larger thermodynamic pressure p = e/3, a peak in the bulk viscosity 
comparable to that indicated by the lattice calculations of Ref. [E7j is sufficient to trigger cavitation. 
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is academic, since the r — > oo solution has P^ negative. In fact, with a power-law divergent (/s the 
bulk shear is large enough to drive P^ < already at rather early times. Cavitation occurs long 
before the asymptotic solution discovered in Ref. [^7[ is reached. 



Following Ref. fl37 |, we have investigated the entropy produced according to the hydrodynamic 
equations that we have solved. For example, comparing the calculation illustrated in the top-left 
panel of Fig. |], in which (/s is as in (34) with ( |3.2|) and (3J3), with the calculation illustrated in 



Fig. [I], in which £ = 0, we find that turning on a bulk viscosity with a peak as in Meyer's lattice 



results [27] increases the entropy by about 20%. This suggests that entropy production due to the 



bulk viscosity is a modest effect, as was found previously in Ref. |S7| for cases in which the peak 
in (/s is not high enough to cause cavitation. In our case, however, this result cannot be trusted 
because there is no way to use our boost-invariant calculation to estimate how much further entropy 
is produced upon cavitation. 

We now consider the effects of increasing r^, which is expected to be large where ( peaks. We 
will not attempt to model a time- varying r^; instead, we ask what are the consequences of increasing 
relative to (|L4|) at all temperatures. If we increase t£ by a factor of 10, keeping AT as in (|3.3|) , 
we find that a sta bie decreases from 0.37 to 0.33. If we increase t£ by a further factor of two, meaning 
that it is 20 times greater than in ( fTiP , a sta bi e increases to 0.49. If we increase by one further 
factor of two, making it 40 times greater than in ( |3.4j ), a sta bie increases to 0.94, comparable to the 
a in ( |3.3| ). In the bottom-right panel of Fig. |4] we illustrate the case with r£ = 20 r^. The bulk 
stress II changes more gradually as a function of time, and as a result P^ does not go as far negative 
as in the top-left panel of Fig. [|, but the change is not dramatic and as a consequence a sta bi e is 
greater, but not by much. 7 This insensitivity to large changes in the value of is one indication 
that second order effects are still small at the time when cavitation occurs. Another indication of 
this is the fact that |n|/(£ + p) is small, for example less than 10% at the time when cavitation 
occurs in the top-left panel of Fig. 4. 8 

We can summarize these results as follows: 



With a and AT as in ( |3.3| ), the peak in the bulk viscosity above T c drives the hydrodynamic 
evolution to negative Pg, indicating cavitation. 

Stable hydrodynamic evolution all the way down to T = T c requires that the peak in (/s 
near T c be less than a threshold that is one quarter to one half as high as the peak found in 



Ref. [27 1; the threshold peak height is fairly insensitive to the width of the peak, for widths 



between one quarter and four times that found in Ref. 

The effects of a peak in (/s near T c can be washed out by increasing the relaxation time 
for the bulk stress II. However, the increase must be by a very large factor. If is larger 
than Tn of fl2.19| ) by a factor of 10 (or 20), results are little affected and the peak in (/s 



7 After the first version of this paper appeared, Song and Heinz found that cavitation does n ot o ccur anywhere 
in their (3 + l)-dimensional calculation if they use = (120 fm)(C/s) [|o), which is larger than (3^) by a factor of 



several hundred in the vicinity of the peak in the bulk viscosity. Our results are consistent with this. 

8 As noted in Section 2.1, once we break conformality by introducing a nonzero bulk viscosity further second order 
terms can arise in both ( [2.12| ) and ( 2.13 ) [Q. See Ref. (^9| for one example. We have confirmed that adding the 



second order terms considered in Ref. ]69| has only negligible effects on our results. 
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must still be reduced by about a factor of three (or two) relative to (3^3) in order to obtain 
stable hydro dynamic evolution all the way down to T = T c . If r£ is 40 times greater than Tq, 
however, the effects of a bulk viscosity peak are washed out sufficiently that P^ remains just 
barely positive even for a peak whose height and width are as in (|3.3|) . 



We have focused on Pg of ( |2.10| ) rather than P± of ( |2.9| ) because the nonzero shear stress $ 
implies that as the bulk stress II becomes increasingly negative, P^ goes negative first, P± only 
later. By the time P± goes negative in the calculation, the calculation has already broken down at 
the time when P^ went negative, triggering cavitation. But, we see in all the panels in Fig. |] that 
$ is quite small by the time P^ goes negative, meaning that P± is already close to zero when Pg 
reaches zero. It will be interesting to see, therefore, which component of the pressure goes negative 
first at which location in the fluid in a calculation that includes transverse expansion. 



4. Implications 

In thinking through the implications of our results, the first possibility to consider is that the peak 
in (/s near T c in QCD with quarks is in fact not so high that cavitation results. As we discussed in 
Section 3.1, the current lattice calculations of the height of the bulk viscosity peak in QCD without 
quarks have various caveats, meaning that the peak in this theory could be smaller (although it 
could just as well be larger) than is indicated by the results of Ref. |2"7fl . Furthermore, there are 



various indications that (/s is somewhat lower in QCD with quarks than in QCD without quarks. 
At very high temperatures where the quark-gluon plasma is weakly coupled, if one compares the 
two theories at a fixed small value of the QCD coupling, say o:qcd = 0.2, one finds that (/s in 
QCD with three flavors of quarks is about 56% of that in quarkless QCD |58| . At these high 
temperatures, (/s is smaller than 10~ 3 in value, so this comparison can give only rough guidance to 
how the height of the peak near T c will change when quarks are introduced, but it does suggest that 
it will decrease. A second argument is simply that the transition in QCD with quarks is a crossover 
whereas that in quarkless QCD is first order, and if adding quarks smooths out the transition then 
it is reasonable to guess that it will also round off the peak in (/s. The magnitude of this effect 
can be guessed by looking at (e — 2>p)/{e + p) which, like (/s, is a dimensionless measure of the 
breaking of conformality. In quarkless QCD, (e — 3p) /(e + p) peaks at a value of 0.85 [BSJ while 



in the equation of state ( |2.16| ) for QCD with quarks, (e — 3p)/(s + p) peaks at 0.53. So, both this 
argument and the comparison to what happens at very high temperatures suggest that the peak in 
(/s in QCD with quarks is (very roughly) about half as high as in (|3.1|) with a as in ( |3.3| ). Taking 
our results at face value, this would put it just above a sta bie, meaning that the peak in (/s would 
trigger cavitation very close to T c . The uncertainties are large and it could certainly be that the 
peak height ends up lower than a s t a bi e) and no cavitation occurs near T c . The previous studies of 
the effects of the peak in the bulk viscosity in Refs. [EJ71 O, [L3| have explored the consequences 



of peaks that are not high enough to cause cavitation. 9 If this is the path that nature chooses, 



9 It is worth noting that in the examples of nonconformal plasmas in which the authors of Refs. j7(| |65| [n], 
have been able to compute (,/s via gauge/gravity duality, a peak in £/s is seen but it is not large enough to cause 
cavitation. For example, in the model of Ref. |}5) the ratio (/r/ is given by 2(-| — c 2 & ), with c s the speed of sound, 
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then hydrodynamic calculations can be followed down to temperatures below T c , and it becomes 
interesting to investigate the possibility of cavitation at a lower temperature, driven by the rising 
shear viscosity of the hadronic phase at low temperatures |23], p4| . 

It is more interesting to consider the possibility that the peak in (/s near T c in QCD with 
quarks is large enough to cause the expanding fluid produced in heavy ion collisions to cavitate 
when it cools through T ~ T c , as our 1 + 1-dimensional calculations indicate. There are a variety 
of aspects of the observed phenomenology of heavy ion collisions that give some support to this 
possibility: 



The possibility that the peak in the bulk viscosity near T c can cause the fluid to fragment, 
and then freezeout, has been considered previously in Refs. JJ9J. These authors suggest that 
data on two particle momentum correlations (the HBT effect) in heavy ion collisions at RHIC 
can be understood if the hadrons in the final state come from such fragments. They have also 
suggested other experimental observables in Ref. [[74]. 



One of the workhorses of heavy ion phenomenology is the statistical hadronization model, 



reviewed in Ref. |75| , which has been used successfully to describe the ratios among the yields 
of many different hadrons using a few parameters including the chemical freezeout temperature 
and baryon chemical potential. One of the conceptual underpinnings of this model, described 



in Ref. |75[ and going back to the original formulation of Hagedorn |76|, is the assumption 
that high energy collisions give rise to multiple clusters — colorless, extended, massive objects 
- which then hadronize statistically (meaning that all hadronic final states consistent with 
conservation laws are equally likely). Most work in this context has focused on the statistical 
hadronization; the dynamics that results in the generation of the clusters in the first place 
has received less attention. This dynamics is no doubt complex, and may be quite different in 
hadron-hadron and heavy ion collisions. Our work suggests that in the case of ultrarelativistic 
heavy ion collisions, in which a hydrodynamic description in terms of an expanding fluid with 
T > T c is appropriate in the early stages of the collision, the clusters required by the statistical 
hadronization model may arise via cavitation, and this cavitation may be triggered by the 
peak in C/s in the vicinity of T c . Perhaps this could be an explanation of why the chemical 
freezeout temperatures extracted via the use of the statistical hadronization model seem to 
be in the vicinity of T c M . 



• In Refs. J77]], Broniowski et al have explained the event-by-event fluctuations in the mean 
transverse momentum of the particles produced in heavy ion collisions at RHIC via the same 
assumption that underlies the statistical hadronization model, namely that hadronization is 
preceded by the material produced in a heavy ion collision falling apart into clusters, each of 
which then yield 6 to 15 hadrons when they hadronize. 

and is therefore everywhere less than 2/3. In contrast, in the example analyzed in Ref. j73| via gauge/gravity duality 
C can be 3> r\ and <^/s rises comparably high to the peak found in Meyer's lattice calculations [ p7| , more than high 
enough to trigger cavitation. At present, these calculations taken together therefore support the existence of a peak 
in £/s but do not provide sufficient guidance regarding its height. 
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• In Ref. f75| , the PHOBOS collaboration provides evidence (from two-particle correlations in 
pseudorapidity and azimuthal angle) that the hadrons in the final state produced in heavy 
ion collisions at RHIC come from clusters that decay into 3 to 6 charged hadrons, meaning 5 
to 9 hadrons in all. 

As we discussed in Section 1, it is pleasing to have a means by which a hydrodynamic calculation 
can predict its own break down. The peak in the bulk viscosity near T c can provide a simple and 
elegant means: if this peak is high enough — as we have quantified — it drives the longitudinal 
pressure to zero at which point the fluid cavitates. The phenomenological evidence in support of 
the notion that hadronization in ultrarelativistic heavy ion collisions is preceded by cavitation, with 
the fluid fragmenting into droplets that play the role of the clusters which have long been employed 
in the statistical hadronization framework, is perhaps not yet overwhelming. But, together with 
our investigation, it certainly seems sufficient to take this possibility seriously. 
There are many avenues open for further investigation: 

• An investigation of the effects of other possible second order terms that can arise in the 
hydrodynamic equations for a nonconformal fluid, as well as third order terms, would be 
desirable. Little is known about the coefficients of such terms. But, although they will have 
quantitative effects, given the insensitivity of our results to large variations in r^, Ai and 
there is no reason to expect qualitative effects. 



A more interesting direction to pursue is to repeat our study using a 3 + 1-dimensional 
hydrodynamic code, describing both longitudinal and transverse expansion. This will make 
the determination of the height of the peak in (/s that is needed in order to trigger cavitation 
when the fluid cools through T c more quantitative. And, having a fluid whose energy density 
varies with tranverse position will raise new questions and open new possibilities. For example, 
cavitation should occur much earlier at the (cooler) edges of the collision region, since T ~ 
T c there earlier [|(| |79 |. Cavitation should start at the edges and move inward, just as 
hadronization has long been understood to do. 



• It is important to investigate whether if freezeout and hadronization are triggered by cavitation 
near T c this modifies the extraction of rj/s via comparison between data and hydrodynamic 
calculations of the anisotropic expansion of the fluid produced in collisions with a nonzero 
impact parameter. The effect of the physics of cavitation near T c on this comparison may 
prove minimal, since the anisotropic flow is generated early in the collision, when the hot fluid 
is still azimuthally anisotropic in shape. This means that the anisotropic flow is generated 
well before cavitation is triggered, when the bulk viscosity is still small compared to the shear 
viscosity. 



From the point of view of understanding the observable, and perhaps observed, phenomenology 
of cavitation, the most important question is the determination of the size distribution of the 
droplets formed when the fluid cavitates. Work in this direction can be found in Refs. [BUI. A 



from first principles determination of the size distribution will be challenging: for one, it will 
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require determining the surface tension associated with the interface between expanding quark- 
gluon plasma with T ~ T c and vacuum. If this surface tension is small, small droplets will be 
favored. The requirement that the droplets must be color singlets will require rearrangement 
of color within ~ Aqq D of the surface as droplets separate during cavitation. Although hard 
to quantify, this can be thought of as a contribution to the surface tension which sets a limit 
on the smallness of the droplets that is of order a few times Aq^ d . It is then interesting to 
note that a spherical droplet with a radius of 1 fm that has the energy density obtained from 
the equation of state ( |2.16|) at T = T c contains about 6 GeV of energy, which is in the right 
ballpark to explain the PHOBOS data J78J mentioned above. This suggests that the surface 
tension is small enough that cavitation yields many small droplets. The picture to have in 
mind is that as the quark-gluon plasma produced in an ultrarelativistic heavy ion collision 
expands and cools through T ~ T c , the fluid falls apart into a mist of hundreds of small 
droplets, each of which later hadronizes as in the statistical hadronization model. Cavitation 
into a mist of small droplets which then become a hadron gas is not likely to have dramatic 
observable consequences 
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